Обнаружение статистически значимых отличий в уровнях экспрессии генов больных раком

Это задание поможет вам лучше разобраться в методах множественной проверки гипотез и позволит применить ваши знания на данных из реального биологического исследования.

В этом задании вы:

вспомните, что такое t-критерий Стьюдента и для чего он применяется сможете применить технику множественной проверки гипотез и увидеть собственными глазами, как она работает на реальных данных почувствуете разницу в результатах применения различных методов поправки на множественную проверку Основные библиотеки и используемые методы:

Библиотека scipy и основные статистические функции:http://docs.scipy.org/doc/scipy/reference/stats.html#statistical-functions

Библиотека statmodels для методов коррекции при множественном сравнении:

http://statsmodels.sourceforge.net/devel/stats.html

Статья, в которой рассматриваются примеры использования statsmodels для множественной проверки гипотез:

http://jpktd.blogspot.ru/2013/04/multiple-testing-p-value-corrections-in.html

Описание используемых данных

Данные для этой задачи взяты из исследования, проведенного в Stanford School of Medicine. В исследовании была предпринята попытка выявить набор генов, которые позволили бы более точно диагностировать возникновение рака груди на самых ранних стадиях.

В эксперименте принимали участие 24 человек, у которых не было рака груди (normal), 25 человек, у которых это заболевание было диагностировано на ранней стадии (early neoplasia), и 23 человека с сильно выраженными симптомами (cancer).

Ученые провели секвенирование биологического материала испытуемых, чтобы понять, какие из этих генов наиболее активны в клетках больных людей.

Секвенирование — это определение степени активности генов в анализируемом образце с помощью подсчёта количества соответствующей каждому гену РНК.

В данных для этого задания вы найдете именно эту количественную меру активности каждого из 15748 генов у каждого из 72 человек, принимавших участие в эксперименте.

Вам нужно будет определить те гены, активность которых у людей в разных стадиях заболевания отличается статистически значимо.

Кроме того, вам нужно будет оценить не только статистическую, но и практическую значимость этих результатов, которая часто используется в подобных исследованиях.

Диагноз человека содержится в столбце под названием "Diagnosis".

Практическая значимость изменения

Цель исследований — найти гены, средняя экспрессия которых отличается не только статистически значимо, но и достаточно сильно. В экспрессионных исследованиях для этого часто используется метрика, которая называется fold change (кратность изменения). Определяется она следующим образом:

Fc(C,T)= T/C, if T > C; -C/T, if C > T.

где C,T — средние значения экспрессии гена в control и treatment группах соответственно. По сути, fold change показывает, во сколько раз отличаются средние двух выборок.

Инструкции к решению задачи

Задание состоит из трёх частей. Если не сказано обратное, то уровень значимости нужно принять равным 0.05.

Часть 1: применение t-критерия Стьюдента

В первой части вам нужно будет применить критерий Стьюдента для проверки гипотезы о равенстве средних в двух независимых выборках. Применить критерий для каждого гена нужно будет дважды:

для групп normal (control) и early neoplasia (treatment) для групп early neoplasia (control) и cancer (treatment) В качестве ответа в этой части задания необходимо указать количество статистически значимых отличий, которые вы нашли с помощью t-критерия Стьюдента, то есть число генов, у которых p-value этого теста оказался меньше, чем уровень значимости.

Часть 2: поправка методом Холма

Для этой части задания вам понадобится модуль multitest из statsmodels.

 1 import statsmodels.stats.multitest as smm В этой части задания нужно будет применить поправку Холма для получившихся двух наборов достигаемых уровней значимости из предыдущей части. Обратите внимание, что поскольку вы будете делать поправку для каждого из двух наборов p-value отдельно, то проблема, связанная с множественной проверкой останется.

Для того, чтобы ее устранить, достаточно воспользоваться поправкой Бонферрони, то есть использовать уровень значимости 0.05 / 2 вместо 0.05 для дальнейшего уточнения значений p-value c помощью метода Холма.

В качестве ответа к этому заданию требуется ввести количество значимых отличий в каждой группе после того, как произведена коррекция Холма-Бонферрони. Причем это число нужно ввести с учетом практической значимости: посчитайте для каждого значимого изменения fold change и выпишите в ответ число таких значимых изменений, абсолютное значение fold change которых больше, чем 1.5.

Обратите внимание, что

применять поправку на множественную проверку нужно ко всем значениям достигаемых уровней значимости, а не только для тех, которые меньше значения уровня доверия. при использовании поправки на уровне значимости 0.025 меняются значения достигаемого уровня значимости, но не меняется значение уровня доверия (то есть для отбора значимых изменений скорректированные значения уровня значимости нужно сравнивать с порогом 0.025, а не 0.05)! Часть 3: поправка методом Бенджамини-Хохберга

Данная часть задания аналогична второй части за исключением того, что нужно будет использовать метод Бенджамини-Хохберга.

Обратите внимание, что методы коррекции, которые контролируют FDR, допускает больше ошибок первого рода и имеют большую мощность, чем методы, контролирующие FWER. Большая мощность означает, что эти методы будут совершать меньше ошибок второго рода (то есть будут лучше улавливать отклонения от H0, когда они есть, и будут чаще отклонять H0, когда отличий нет).

В качестве ответа к этому заданию требуется ввести количество значимых отличий в каждой группе после того, как произведена коррекция Бенджамини-Хохберга, причем так же, как и во второй части, считать только такие отличия, у которых abs(fold change) > 1.5.


In [23]:
import numpy as np
import pandas as pd
import statsmodels
import scipy as sc
import statsmodels.stats.multitest as smm

In [3]:
data = pd.read_csv('gene_high_throughput_sequencing.csv')

In [4]:
data.head()


Out[4]:
Patient_id Diagnosis LOC643837 LOC100130417 SAMD11 NOC2L KLHL17 PLEKHN1 C1orf170 HES4 ... CLIC2 RPS4Y1 ZFY PRKY USP9Y DDX3Y CD24 CYorf15B KDM5D EIF1AY
0 STT5425_Breast_001_normal normal 1.257614 2.408148 13.368622 9.494779 20.880435 12.722017 9.494779 54.349694 ... 4.761250 1.257614 1.257614 1.257614 1.257614 1.257614 23.268694 1.257614 1.257614 1.257614
1 STT5427_Breast_023_normal normal 4.567931 16.602734 42.477752 25.562376 23.221137 11.622386 14.330573 72.445474 ... 6.871902 1.815112 1.815112 1.815112 1.815112 1.815112 10.427023 1.815112 1.815112 1.815112
2 STT5430_Breast_002_normal normal 2.077597 3.978294 12.863214 13.728915 14.543176 14.141907 6.232790 57.011005 ... 7.096343 2.077597 2.077597 2.077597 2.077597 2.077597 22.344226 2.077597 2.077597 2.077597
3 STT5439_Breast_003_normal normal 2.066576 8.520713 14.466035 7.823932 8.520713 2.066576 10.870009 53.292034 ... 5.200770 2.066576 2.066576 2.066576 2.066576 2.066576 49.295538 2.066576 2.066576 2.066576
4 STT5441_Breast_004_normal normal 2.613616 3.434965 12.682222 10.543189 26.688686 12.484822 1.364917 67.140393 ... 11.227770 1.364917 1.364917 1.364917 1.364917 1.364917 23.627911 1.364917 1.364917 1.364917

5 rows × 15750 columns


In [6]:
data.shape


Out[6]:
(72, 15750)

In [8]:
data['Diagnosis'].value_counts()


Out[8]:
early neoplasia    25
normal             24
cancer             23
Name: Diagnosis, dtype: int64

In [9]:
normal = data[data['Diagnosis'] == 'normal']

In [11]:
neoplasia = data[data['Diagnosis'] == 'early neoplasia']

In [13]:
cancer = data[data['Diagnosis'] == 'cancer']

In [33]:
def get_fc(c, t):
    if t >= c:
        return t / c
    else:
        return -c / t

In [34]:
n1, n2 = 0, 0
p1 = np.zeros(data.shape[1] - 2)
p2 = np.zeros(data.shape[1] - 2)
fc1 = np.zeros(data.shape[1] - 2)
fc2 = np.zeros(data.shape[1] - 2)
for i in xrange(2, data.shape[1]):
    gene_normal = normal.iloc[:, i]
    gene_neoplasia = neoplasia.iloc[:, i]
    gene_cancer = cancer.iloc[:, i]

    p1[i-2] = sc.stats.ttest_ind(gene_normal, gene_neoplasia, equal_var = False).pvalue
    p2[i-2] = sc.stats.ttest_ind(gene_neoplasia, gene_cancer, equal_var = False).pvalue
    
    fc1[i-2] = get_fc(gene_normal.mean(), gene_neoplasia.mean())
    fc2[i-2] = get_fc(gene_neoplasia.mean(), gene_cancer.mean())
    
n1 = len(p1[p1[:] < 0.05])
n2 = len(p2[p2[:] < 0.05])
    
print n1, n2
with open('task1.1.txt', 'w') as f:
    f.write(str(n1))
with open('task1.2.txt', 'w') as f:
    f.write(str(n2))


1575 3490

In [39]:
len(fc1[fc1[:] > 1.5]), len(fc2[fc2[:] > 1.5])


Out[39]:
(156, 451)

In [43]:
n1, n2 = 0, 0
reject, p1_corr, a1, a2 = smm.multipletests(p1, alpha = 0.025, method = 'holm') 
reject, p2_corr, a1, a2 = smm.multipletests(p2, alpha = 0.025, method = 'holm')

for i in xrange(len(p1_corr)):
    if p1_corr[i] < 0.025 and np.abs(fc1[i]) > 1.5:
        n1 += 1
    if p2_corr[i] < 0.025 and np.abs(fc2[i]) > 1.5:
        n2 += 1
        
print n1, n2
with open('task2.1.txt', 'w') as f:
    f.write(str(n1))
with open('task2.2.txt', 'w') as f:
    f.write(str(n2))


2 77

In [44]:
n1, n2 = 0, 0
reject, p1_corr, a1, a2 = smm.multipletests(p1, alpha = 0.025, method = 'fdr_bh') 
reject, p2_corr, a1, a2 = smm.multipletests(p2, alpha = 0.025, method = 'fdr_bh')

for i in xrange(len(p1_corr)):
    if p1_corr[i] < 0.025 and np.abs(fc1[i]) > 1.5:
        n1 += 1
    if p2_corr[i] < 0.025 and np.abs(fc2[i]) > 1.5:
        n2 += 1
        
print n1, n2
with open('task3.1.txt', 'w') as f:
    f.write(str(n1))
with open('task3.2.txt', 'w') as f:
    f.write(str(n2))


4 524

In [ ]: